Bifurcations in the optimal elastic foundation for a buckling column 
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Abstract 

We investigate the buckling under compression of a slender beam with a distributed lateral elastic support, for which 
there is an associated cost. For a given cost, we study the optimal choice of support to protect against Euler buckling. 
We show that with only weak lateral support, the optimum distribution is a delta-function at the centre of the beam. 
When more support is allowed, we find numerically that the optimal distribution undergoes a series of bifurcations. 
We obtain analytical expressions for the buckling load around the first bifurcation point and corresponding expansions 
for the optimal position of support. Our theoretical predictions, including the critical exponent of the bifurcation, are 
confirmed by computer simulations. 



■ 1. Introduction 

^ ■ n 

Buckling is a common mode of mechanical failure 
q and its prevention is key to anysuccessful engineering de- 
^ sign. As early as 1759, Euler [2[ gave an elegant descrip- 
tion of the buckling of a simple beam, from which the 
" 55 so-called Euler buckling limit was derived. Works which 
>~>cite the goal of obtaining structures of least weight stable 
i-C iagainst buckling can be found throughout the literature 
l^J'S EL Hi S 0' an d much understanding has been gained 
on optimal structural design [H d, 0] . Designs of ever in- 
creasing complexity have been analysed and recent work 
suggested that the optimal design of non-axisymmetric 
columns may involve fractal geometries 13, 11|- With the 
development of powerful computers, more and more com- 
plicated structures can be designed with optimised me- 
chanical efficiency. However, understanding and prevent- 
ing buckling remains as relevant as ever. 

In this paper, we consider a simple uniform elastic 
beam, freely hinged at its ends and subjected to a compres- 
sive force and therefore vulnerable to buckling. However, 
in contrast to Euler's original problem, we specify that the 
beam is stabilized by restoring forces, perpendicular to its 
length, which are provided by an elastic foundation (as 
illustrated in figure [1]). This represents a simple and prac- 
tical method of protecting against buckling instabilities. 

In the simplest case figure l(a)| we can imagine this 
elastic foundation as a finite collection of linear springs at 
points along the beam. Each has a spring constant, and so 
provides a restoring force at this point, proportional to the 
lateral deflection of the beam. More generally, the elastic 
foundation could be distributed as a continuous function 
along the length of the beam, rather than being concen- 
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is a spring constant per unit length, which may vary along 
the beam. 

We are interested in optimising this elastic support, 
and so we need to specify a cost function for it. This 
we take to be the sum of the spring constants (if there 
are a discrete collection of springs) or the integral over 
the spring constant per unit length along the beam (if the 
elastic foundation is continuous) . By choosing the optimal 
distribution of these spring constants, we wish to find the 
minimum cost of elastic support which will protect against 
buckling under a given compressive load (or equivalently, 
the distribution of an elastic support of fixed cost which 
will support the maximum force). 

The optimal position of one or two deformable or in- 
finitely stiff supports have been studied in the literature 
(see for example Rcf. 12] and references therein), and gen- 
eral numerical approaches established for larger numbers 
of supports 12[. However, in the present paper, we con- 



trated into discrete springs (figure 1(b) I. In this case, there 



sider the general case where any distribution of support is 
in principle permitted. 

A perturbation analysis shows that in the limit of weak 
support strength, the optimal elastic foundation is a con- 
centrated delta-function at the centre of the beam, but 
when stronger supports are permitted, we show that the 
optimal solution has an upper bound on the proportion of 
the beam that remains unsupported. In this sense, the op- 
timum distribution becomes more uniform for higher val- 
ues of support strength. To tackle the problem in more 
detail, we develop a transfer matrix description for the 
supported beam, and we find numerically that the opti- 
mal supports undergo a series of bifurcations, reminiscent 
of those encountered in iterated maps. However, we are 
only able to proceed a limited distance in the parameter 
space and we are unable to explore for more complex be- 
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haviour (for example, any possible signature of chaos |13l|). 

We obtain analytic expressions for the buckling load 
in the vicinity of the first bifurcation point and a corre- 
sponding series expansion for the optimal placement of 
elastic support. Following this optimization we show that 
a mathematical analogy between the behaviour exhibited 
in this problem and that found in Landau theory of second 
order phase transitions [l3| exists. However, the analogue 
of free energy is non-analytic, while in Landau theory it 
is a smooth function of the order parameter and the con- 
trol variable. Our results, including critical exponents are 
confirmed by computer simulations, and should provide a 
basis for future analysis on higher order bifurcations. 

2. Theory 

A slender beam of length L, hinged at its ends, under 
a compressive force F, is governed by the Euler-Bernoulli 
beam equation [![: 



ay d v 
EI-£ + F + q (x) = o, 



dx 4 



dx 



(1) 



where E is the Young modulus of the beam, / is the second 
moment of its cross sectional area about the neutral plane, 
y is the lateral deflection, x the distance along the beam 
and q(x) is the lateral force applied per unit length of 
beam. The beam is freely hinged at its end points and 
therefore the deflection satisfies y = y" = at x — and 
L. 

If the lateral force is supplied by an elastic foundation, 
which provides a restoring force proportional to the lat- 
eral deflection, then through rescaling we introduce the 
following non-dimensional variables x = ttx/L, y = ny/L, 
f = FL 2 /{EIir 2 ) and p = qL 4 / '(EIir 4 y). Eq. © becomes 



d 4 y d 2 y 

d^ + f dx^ +P{x)v 



for sce(0,7r), (2) 



where y(Q) = y(n) — y"(0) — y"(n) = and p(x) repre- 
sents the strength of the lateral support (for example the 
number of springs per unit length) at position x. 

We are always interested in the minimum value / m i n 
of / that leads to buckling [in other words, the smallest 
eigenvalue of Eq. ([2"])]. For the case of no support (p = 
0), the possible solutions to Eq. @ are / £ Z + , and so 
buckling first occurs when / = 1. 

Lateral support improves the stability (increasing the 
minimum value of the applied force / at which buckling 
first occurs), but we imagine that this reinforcement also 
has a cost. In particular, for a given value of 



p{x) dx, 



(3) 



we seek the optimal function p(x) which maximises the 
minimum buckling force / m i n . 

The simplest choice we can imagine is that p takes the 
uniform value m/ir, so that the form of deflection is y{x) oc 
sin kx, for some integer fc, which represents a wavenumber. 
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Figure 1: Schematic of a slender beam with elastic support, loaded 
under compression force F. (a) shows the case where the lateral 
restoring force per unit length along the beam q(x) is provided by 
linear springs of spring constant at discrete points so that 

= 'Ylt—l kiS{x — Xi)y/L. (b) shows schematically the case where 
there is a continuous support: the lengths of the arrows indicate the 
local spring constant. 



This leads immediately to the result that in this case 



/ min = min 



m 
irk 2 



(4) 



Eq. ((4]) has a physical interpretation: the first term 
comes from the free buckling of the column which is most 
unstable to buckling on the longest allowed length scales 
(i.e. the smallest values of k), as demonstrated by Eu- 
ler. The second term represents the support provided by 
the elastic foundation, which provides the least support at 
the shortest length scales (largest values of fc). The bal- 
ance between these two terms means that asm-> oo, the 
uniformly supported column buckles on a length scale of 
approximately 



^cff ~ (Tr/m) 1 / 4 as m 



oo, 



and can support a load 



/ uni ~ 



(5) 



(6) 



Now, although a uniform elastic support is easy to anal- 
yse, it is clear that this is not always optimal. Consider 
the case where m is very small, so that p provides a small 
correction in Eq. @. In this case, the eigenvalues remain 
well-separated, and we can treat the equation perturba- 
tively: let 



y = y sin x + yi(x) and / = 1 + fx, 



(7) 



then from Eq. ©, if we multiply through by sinx (the 
lowest unperturbed eigenfunction) and integrate, we have 
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to leading order: 

~d 4 yi , d 2 yi 



where T n is given by 









sin a; 







dx 4 dx 2 



?/o sin 2 x[p- fx] 



dx = 0. 



(8) 

Repeated integrations by parts with the boundary condi- 
tions y'( — at x — 0,7r establishes the self-adjointness of 
the original operator, and we arrive at 



fi 



2 f n 

— / p(x) sin x dx. 
71 Jo 



(9) 



We therefore see that in the limit m — > 0, the optimal 
elastic support is p{x) = mS(x — 7r/2), and for this case, 
/ min = l + (2m/7r) + 0(m 2 ). 

The requirement for optimal support has therefore con- 
centrated the elastic foundation into a single point, leaving 
the remainder of the beam unsupported. 

3. Transfer Matrix formulation 

In order to proceed to higher values of m in the opti- 
mization problem, we assume that there are N — 1 discrete 
supports at the positions {x n }, with corresponding set of 
scaled spring constants {/?„}, adding up to the total m: 



p(x) 
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(15) 



Ax n 

S n 
K n 



pl/2 



sin[/ /2 (x n - x n _D 
cos[f 1/2 (x n - X n -i)]. 



At the two end-points at x = 0, 7r, we have the bound- 
ary conditions that y and y" vanish, which leads to the 
following four conditions 



B = D 
An-iSn + Bn-iKn 
Cn-i(xn - xn-i) + -Djv-i 

If we now define a matrix 

R = ■ 
then Eqs. (|16I18|) lead to 

M 



= 
= 
= 0. 

.T 2 T! 



A 
C 



= 



These discrete supports divide the beam into N (not nec- 
essarily equal) segments, and for convenience in later cal- 
culations, we also define the end points as xq = and 

Xn = 7T. 

For each segment of the beam given by x n < x < x n +i, 
the Euler-Bernoulli equation @ can be solved in the form 

y(x) = A n sin[/ 1/2 (a; - x n )] + B n cos[/ 1/2 (x - x n )] 

+C n (x -x n )+D n . (12) 

If we integrate Eq. over a small interval around x n , we 
find that, 



where 
M = 



RuSn + R21KN 
(Ax N )R 3 i + Rai 



Rl3$N + R23KN 
(Ax N )R 33 + i? 43 



(16) 
(17) 
(18) 



(19) 



(20) 



(21) 



y ( x n) = y ( x n ) . v = v (x n ) , 

y"(xi)=y"(x-), 
y'"(x+)-y"'(x-)+p n y(x n ) = 0, (13) 

where 2+ and x~ are values infinitesimally greater and less 
than than x n respectively. Defining v„ = (A n , B n , C n , D n ) T ', 
these continuity constraints on the piecewise solution of 
Eq. (|12[) can be captured in a transfer matrix 



For the beam to buckle, there needs to be non-zero solu- 
tions for Ao and/or Cq. Therefore, the determinant of M, 
which is a function of /, must go to zero. The smallest 
/) /min, at which det(M) = 0, gives the maximum com- 
pression tolerated by the beam and its support. The task, 
thus, is to find the set of {/3 n } and {x n } which maximise 

/min- 

4. Equally spaced, equal springs 

Any definite choice of p[x) provides a lower bound on 
the maximum achievable value of f mm , so before discussing 
the full numerical optimization results on p, we consider 
here a simple choice of p which illuminates the physics. 

Suppose that p consists of equally spaced, equally strong 
delta-functions: 



Pn(x) 



N-l 



^ N- 

n=l 



1 



8{x — ir/n). 



(22) 



V„-l, 



(14) 



The value of / m i n can be found by a straight-forward calcu- 
lation for each value of m, using the transfer matrix formu- 
lation above. The results are plotted in figure^ and we see 
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Figure 2: Value of / m j n for p constant (dashed line), and for equally 
spaced, equally strong delta functions (N is the number of intervals, 
so N — 1 is the number of delta-functions). 
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Figure 4: Showing the critical exponents for the first and second 
bifurcation in the restricted problem of /3 n = m/(N — 1). 
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Figure 3: Results of the restricted optimization, obtaining the set 
{x n } with constant n = m/(N — 1). 



that in general, it is better to concentrate the elastic sup- 
port into discrete delta functions, rather than having a uni- 
form elastic support. However, it is important to choose 
the appropriate number of delta functions: if the number 
is too few, then there will always be a buckling mode with 
/ = TV 2 which threads through the comb of delta functions 
without displacing them. However, apart from this con- 
straint, it appears to be advantageous to choose a smaller 
value of N; in other words, to concentrate the support. 

5. Numerical optimization of the support 

Before we look at the general optimization problem 
where we will seek the optimal set of {x n } and {/3 n } for a 
given cost, we investigate a simplified problem to give us 
further insight into the nature of the problem. We set 

777 

and then find the set {x n } which maximises f m in- The 
results obtained from an exhaustive search are shown in 



figure [3j where we find two bifurcation points in the range 
< m < 40. The critical exponent of each has been ob- 
tained through simulation as, 

a x = 0.5 ±0.01, (24) 
a 2 = 0.49 ±0.03, (25) 

for the first and second bifurcation respectively. Figure 2] 
shows the data from which the exponents are taken, where 
values of mo and xq used are, 

m = 5.09, 26.99 (26) 
x = 0.5, 0.281 (27) 

for the first and second bifurcation respectively. The value 
of xo for the lower branching event at m — 26.99 is related 
to the upper branch by symmetry about the midpoint of 
the beam. As discussed previously, the optimal solution 
must split further at higher values of m. We hypothesize 
that within this restricted problem these splits will take 
the form of bifurcations similar in nature to those found 
here. 

Now we turn to the full optimization problem, where 
the values {ft} as well as the positions {xi} of the sup- 
ports may vary. Using the transfer matrix formulation, we 
seek the optimal elastic support consisting of delta func- 
tions. Figure [5] shows the best solutions, found from an 
exhaustive search of four delta functions (N — 5), up to 
m = 50. We see in figure 4 that there are two bifurca- 
tion events, and one coalescence of the branches. Because 
the optimal solution cannot contain long intervals with no 
support (see section [7] below), we expect that if continued 
to larger values of m and N, a series of further bifurcation 
events would lead to a complex behaviour which would 
eventually fill the interval with closely spaced delta func- 
tions as m — > oo. 

6. First branch point 

Numerical results (figure [5]) indicate that although a 
single delta function at x = ir/2 is the optimal form for 
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Figure 5: Value of f m i n for the optimal form of p(x) and also for 
comparison p constant, and for equally spaced, equally strong delta 
functions. 
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Figure 6: Position of optimal springs as a function of m. The area of 
each circle is proportional to the strength ft of the relevant support, 
with the total area of all the circles at each value of m chosen to be 
a constant, independent of m. 
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Figure 7: Three dimensional plot of / m i n as a function of the position 
parameter £ and fi = m — (I6/71-). 



p in the limit m — > 0, at some point the optimal support 
bifurcates. 

It is clear that this first bifurcation must happen at 
/ = 4, since this represents the excitation of the first anti- 
symmetric buckling mode in the unsupported beam, and 
the delta function at x = tt/2 provides no support against 
this mode. Although the value of / at this first branch 
point is clear, neither the value of m at which it occurs, 
nor the nature of the bifurcation are immediately obvious. 

In order to clarify the behaviour at this first branch 
point, we perform a perturbation expansion: Let us sup- 
pose that N = 3 and 



. . m . ( 7T \ 771 / 7T 



(28) 



where £ and — £ are clearly equivalent, and we will quote 
only the positive value later. Thus {xq, xi, X2, X3} are 
given by {0, tt/2 - £, tt/2 + £, tt} and #,.=#, = m/2. 

We wish to evaluate the matrix M in Eq. (|2"Tj) and seek 
the smallest / giving a zero determinant. On performing 
a series expansion of the determinant for / near 4, we find 
that the critical value of m is I6/71". Furthermore, if we 
define small quantities \x and £ through 



m 



16 16 , 

1- U = h/i£ 

7T 7T 



(29) 

Z = e\e\ (30) 
where e< 1 and £' and \j! are order 1 quantities and 

/ = (31) 

then we can perform a series expansion of det(M) in the 
neighbourhood of e = 0, to obtain term by term a series 
expansion for /. We find that there are two solutions, / + 
and /_, which correspond to functions y(x) symmetric and 
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anti-symmetric about x = tt/2 respectively: 

7T 2 2 7T 3 (6-^ 2 ) 



/-t 
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4 + — u zr 
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124416 
7r 4 (27r 2 - 21 
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11943936 



7T 5 (315 - 15tt 2 - n 4 
4299816960 



V + 0( M 6 ) 



-If I [o + o( M 5 )] 



+e 2 
He 3 1 



27T 7T 2 2 7T 3 (3 + 7T 2 ) 

— u h u 2 i '-u, 6 + O(ir) 

9 M 72^ 93312 p VM ; 




o.i 



128 40 
IhT ~ 27 M 



7T 15 - 7T 2 



486 



V + 0( M 3 ) 



+e 4 [o + o(/i 2 )] + ie 5 



1024 

1357T 



(32) 



/_=4+|C|[0 + O(/i 5 )]+e 2 
(128tt 2 + 576) (8 
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+|C 3 |[0 + O(m 3 )] 
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Figure 8: Curve shows the locus of optimal values for £ near the first 
bifurcation point. This divides the £ — n plane into three regions, in 
which / m i n is given by either Eq. (1326 or (1331) as indicated. 



deformation of the beam, then the energy of our system is 
given [l| by 



1 

U =2 



d 2 z\ fdz 
dx 2 J \dx 



p(x)z* 



dx. 



(35) 



3tt 3 

-\e\ 



(33) 



The final value for / m i n in this neighbourhood is then 
/min = min (/+,/-)■ 

The results are plotted in figure [3 and we see that 
the behaviour of / m ; n around the bifurcation point is not 
analytic, since the transition between the two branches 
/+ and /_ leads to a discontinuity in the derivatives of 
/min- The maximal value of / m in (i-e the optimum we are 
seeking), occurs for £ = when /i < 0, and along the locus 
/+ = /_ when /i > 0. 

From Eqs. (|32p and (|33p . this leads to the optimal 
value of £ being 



Any deformation z{x) which results in U[z(x)] < means 
that the beam will be energetically allowed to buckle under 
this deflection. Furthermore, the associated value of / 
which just destabilises the system against this deformation 
cannot be smaller than the lowest buckling mode / m in- 
Consider therefore a particular choice for z[x), namely 



z(x) 





sin 2 




TZ(X — X\) 

X2—X1 



x 6 (0, X\) 

X e [X1,X2) 
X G [X 2 , 1) 



(36) 



which vanishes everywhere except on the interval = 
(x±,X2), which is of length A = X2 — x±. Then Eq. (|3"5")l . 
together with the observation above about / m in leads to 



2tt 2 



opt 



g^/ 2 _^ + 0(M 3/ 2) ifM > Q 

if n < 



fm.w[p(x)} <^- + ^ 



p(x) sin 



7r(x — xx) 



A 



dx. (37) 



(34) 



This is shown in figure together with the regions of the 
(j, — £ plane in which /+ and /_ apply. 

7. Limit of large support stiffness 

The results of our numerical optimisation suggests that 
the optimum support continues to take the form of a dis- 
crete set of delta-functions. Here we investigate the possi- 
ble form of the optimal support in the limit of large m. 

As m increases, the optimal distribution function p opt 
must become more evenly distributed over the interval. To 
see in what sense this is true, we note that the eigenvalue 
problem for buckling modes given by Eq. @ can be de- 
rived from an energy approach: Suppose that z(x) is any 



IT 2 , 

n Jn 

Trivially, we note from the definition of p opt , that 

Vp : fmin\fi(x)] < fmm[Papt(x)], (38) 

so that from Eqs. ([5]), (13"T|) and (|3"5|) . we finally arrive at a 
condition for how evenly distributed p opt must be for large 
in: 

27r 3/2 m l/2 2 ^4 



p pt(x) sin 



n(x — x%) 
A 



da; > 



A 



A 3 ' 

(39) 

A simple corollary of Eq. (|39p is that if p opt is zero on any 
interval fl of length A, then it must be the case that 

A < tt 5 / 4 ™- 1 / 4 . (40) 

The scaling of this length with m is the same as the effec- 
tive buckling length of a uniformly supported beam dis- 
cussed earlier. 
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8. Discussion 



9. Acknowledgements 



The optimal elastic support for our column appears to 
display complex behaviour: at small values of m the sup- 
port is a single delta function, and even at large values 
of m, it appears to be advantageous for p(x) to be con- 
centrated into discrete delta-functions rather than to be a 
smooth distribution. 

Furthermore, the manner in which the system moves 
from a single to multiple delta functions is not trivial, 
and appears to be through bifurcation events. In the 
full optimization problem we find that the first bifurcation 
event occurs with critical exponent of one half. Inverting 
Eq. (J3U) and substituting it into either Eq. ([22]) or ^/ we 
find that, 



fn 



32 c2 

w 2 Sopt _ 3^4 



l 4 + f^ik^ 

while to leading order, 



if n < 0. 



(41) 



pi 



= 1 8V5 

lo 



3/ y /2 



if /i > 
if fi < 0. 



(42) 



In this form, the mathematical similarities to Landau the- 
ory of second order phase transitions become apparent, 
with £ opt playing the role of the order parameter, \i the 
reduced temperature and —f m in the free energy to be min- 
imized. 

However, there is an important difference. In Landau 
theory of second order phase transitions, the free energy 
.Flan is assumed to be a power series expansion in the order 
parameter tp with leading odd terms missing: 



Fan = F + a 2 ^ 2 + a 4 ^ 4 + 



(43) 



where a 2 oc (T—T c ), the reduced temperature. In our case, 
the buckling force / has to be first optimised for even and 
odd buckling. Thus — / m ; n (which is the analogue of -Flan) 
is a minimum over two intersecting surfaces (figure [7]) and 
so non-analytic at the point of bifurcation. 

Nevertheless, the mathematical form of the solution in 
Eq. (|4"2")1 is the same, including the critical exponent. Fur- 
thermore, our numerical results show that, for the equal 
support case, the critical exponent a is preserved for the 
next bifurcation, suggesting that the nature of subsequent 
bifurcations will also remain the same. 

The details of the behaviour for larger values of m is 
as yet unclear: we speculate that there will be a cascade 
of bifurcations, as seen in the limit set of certain iterated 
maps [l5f ; it remains an open question whether there is an 
accumulation point leading to potential chaotic behaviour. 

Further investigation of this regime may shed light on 
structural characteristics required to protect more complex 
engineering structures against buckling instabilities. 
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